To explore this we used the following parameters space :
constant_env=getAlllSummaries("unifPop50k")
I was not keeping track of the outpute rate so the convertion in real generation has to be done by hand:
constant_env$outputrate=10
constant_env=updateScale(constant_env)
Check the trajectories of the simulation for the population size for different mutation rate \(\mu\)(from left to right: \(\mu = 0, \mu = 0.001,\mu=0.01\)), different selective pressure \(\sigma_s\) (from top to bottomright: \(\sigma_s = 0.1, \sigma_smu = 0.2,\sigma_s=0.4,\sigma_s=1000\))and different carrying capacity \(K\)
#pdf("trajectory_N_const.pdf",width=12,height=16)
plotAllTrajVar(constant_env,m=0.05,E=0,obs="N")
#dev.off()
#pdf("trajectory_var_x_const.pdf",width=12,height=16)
plotAllTrajVar(constant_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))
#dev.off()
Check the variance at the equilibrium:
#pdf("equilibrium_var_x_const.pdf",width=12,height=16)
plotAllVariableSummaries(constant_env,E=0,estimate=eq2833b) #we use only simulations without noise
#dev.off()
we use different \(delta\) to generate environment with noise of growing variance
d=c(0.00,0.01,0.10,0.20,0.40,1.00,10.00)
cols=colorRampPalette(c("black","grey"))(length(d))
names(cols)=d
omegas=1:3
plot(1,1,xlim=c(1,500),ylim=c(-10,10),type="n",xlab="t",ylab=expression(theta[t]))
for ( i in rev(d))
lines(1:500,environment(500,omega=0,delta=i),col=cols[as.character(i)])
legend("topleft",legend=paste0("v=",d),col=cols,lty=1,bg="white")
To explore this we used the following parameters space:
v=getAlllSummaries("growingDelta",exclude="outputrate")
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
noisy_env = rbind(v,constant_env[,colnames(v)])
noisy_env = noisy_env[ noisy_env$mu >0,]
Check the trajectory off the effective population size:
#pdf("trajectory_N_noise.pdf",width=12,height=16)
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="N")
#dev.off()
Check the trajectory off the mean variance size:
#pdf("trajectory_var_x_noise.pdf",width=12,height=16)
plotAllTrajVar(noisy_env,m=0.05,E=0,obs="var_x",ylim=c(0,.005))
#dev.off()
Check and compare the final variance :
#pdf("equilibrium_var_x_noise.pdf",width=12,height=16)
plotAllVariableSummaries(noisy_env,E=0,estimate=eq2833b)
#dev.off()
Some of the data seems to be missing because when \(\delta>\sigma\) we have extinctions, has shown in the graph:
#pdf
sigmas=unique(noisy_env$sigma)
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env[ noisy_env$sigma==0.4 & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == 0))
for(s in sigmas){
tmpc=noisy_env[ noisy_env$sigma==s & noisy_env$mu == 0.01 & noisy_env$m == 0.2 ,]
Ks=sort(unique(tmpc$K))
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$delta,mean)
lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("right",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
#dev.off()
We introduce autocorrelation in the environment. To introduce autocorrelation we increase slightly \(omega\)
par(mfrow=c(1,3))
plot(1:500,environment(500,omega=1,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=2,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
plot(1:500,environment(500,omega=3,delta=.1),type="l",xlab="t",ylab=expression(theta[t]))
For this we re-run previous experiment but with \(\omega \in {1,2}\)
omega1=getAlllSummaries(c("omega1growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega1=updateScale(omega1)
omega1$omega=1
omega2=getAlllSummaries(c("omega2growingDelta"))
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
## Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
omega2=updateScale(omega2)
omega2$omega=2
noisy_env_omega = rbind(v,constant_env[,colnames(v)],omega1[,colnames(v)],omega2[,colnames(v)])
Let’s look at the extinction again for the different omegas
par(mfrow=c(1,3))
omegas=unique(noisy_env_omega$omega)
sigmas=sort(unique(noisy_env_omega$sigma))
deltas=sort(unique(noisy_env_omega$delta))
Ks=sort(unique(noisy_env_omega$K))
for(o in omegas){
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
tmpc=noisy_env_omega[ noisy_env_omega$sigma==0.4 & noisy_env_omega$mu == 0.01 & noisy_env_omega$omega == o & noisy_env_omega$m == 0.2 ,]
xrange=range(noisy_env_omega$delta[noisy_env_omega$omega>0])
plot(tmpc$delta,tmpc$extinction,log="xy",type="n",ylab="extinction time",xlab=expression(delta),main=bquote(mu == 0.01 ~ omega == .(o)),xlim=xrange)
for(s in sigmas){
tmpc=noisy_env_omega[ noisy_env_omega$sigma==s & noisy_env_omega$mu == 0.01 & noisy_env_omega$m == 0.2 & noisy_env_omega$omega == o,]
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$delta,mean)
lines(sort(unique(t$delta)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("bottomleft",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
}
We can thus reproduce he trajectories and variance equilibrium using simulation with \(\omega=1\) et \(\omega=2\)
Check the trajectory off the effective population size:
#pdf("trajectory_N_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="N")
#dev.off()
#pdf("trajectory_N_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="N")
#dev.off()
Check the trajectory off the mean variance size:
#pdf("trajectory_var_x_omega1.pdf",width=12,height=16)
plotAllTrajVar(omega1,m=0.05,E=0,obs="var_x",ylim=c(0,.002))
#dev.off()
#pdf("trajectory_var_x_omega2.pdf",width=12,height=16)
plotAllTrajVar(omega2,m=0.05,E=0,obs="var_x",ylim=c(0,.002))
#dev.off()
And check the variance at the end of the equilibrium:
constant_env$omega=0
omega1_wdz=rbind(omega1,constant_env[constant_env$mu>0,])
omega2_wdz=rbind(omega2,constant_env[constant_env$mu>0,])
omega1_wdz=omega1_wdz[omega1_wdz$sigma<10,]
omega2_wdz=omega2_wdz[omega2_wdz$sigma<10,]
#pdf("equilibrium_var_x_omega1.pdf",width=12,height=16)
plotAllVariableSummaries(omega1_wdz,E=0,estimate=eq2833b) #we use only simulations without noise
#dev.off()
#pdf("equilibrium_var_x_omega2.pdf",width=12,height=16)
plotAllVariableSummaries(omega2_wdz,E=0,estimate=eq2833b) #we use only simulations without noise
#dev.off()
par(mfrow=c(3,3))
deltas=c(0,0.1,1)
vts=c(0.001,0.002,0.003)
cols=colorRampPalette(c("darkgreen","green"))(length(vts))
names(cols)=vts
omegas=1:3
for(d in deltas){
for(o in omegas){
plot(1,1,xlim=c(1,500),ylim=c(-(d+d*.1),1.1*500*max(vts)+d),type="n",xlab="t",ylab=expression(theta[t]),main=bquote(delta == .(d)~omega == .(o)))
for ( vt in vts)
lines(1:500,environment(500,omega=o,delta=d,vt=vt),col=cols[as.character(vt)])
legend("topleft",legend=paste0("v=",vts),col=cols,lty=1)
}
}
moving_theta=getAlllSummaries("movingTheta")
moving_theta=updateScale(moving_theta)
#pdf("extinctionsMovingTheta.pdf",width=17,height=4)
par(mfrow=c(1,3))
rates=sort(unique(moving_theta$vt))
sigmas=sort(unique(moving_theta$sigma))
deltas=sort(unique(moving_theta$delta))
Ks=sort(unique(moving_theta$K))
cols=colorRampPalette(c("darkgreen","yellow"))(length(sigmas))
names(cols)=sigmas
for(d in deltas){
tmpcs=moving_theta[ moving_theta$mu == 0.001 & moving_theta$delta == d & moving_theta$m == 0.1 ,]
#xrange=range(moving_theta$delta)
plot(tmpcs$vt,tmpcs$extinction,log="y",type="n",ylab="extinction time",xlab=expression(v),main=bquote(delta == .(d)~mu == 0.01 ~ omega == .(o)))
for(s in sigmas){
tmpc=tmpcs[ tmpcs$sigma==s ,]
for(k in 1:length(Ks)){
t=tmpc[tmpc$K==Ks[k],]
r=tapply(t$extinction,t$vt,mean)
lines(sort(unique(t$vt)),r,type="b",pch=k,col=cols[as.character(s)])
}
}
legend("bottomleft",
legend=c(paste0("K=",Ks),sapply(sigmas,function(s)as.expression(bquote(sigma==.(s))))),
col=c(rep(1,length(Ks)),cols),
lty=c(rep(NA,length(Ks)),rep(1,length(sigmas))),
pch=c(seq_along(Ks),rep(NA,length(sigmas)))
)
}
#dev.off()
Check the trajectory of the effective population size:
for(t in seq_along(rates)){
#pdf(paste0("trajectory_N_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="N")
#dev.off()
}
Check the trajectory of the mean variance size:
for(t in seq_along(rates)){
#pdf(paste0("trajectory_var_x_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="var_x",ylim=c(0,.005))
#dev.off()
}
Check the trajectory of the distance to optimum:
for(t in seq_along(rates)){
#pdf(paste0("trajectory_dist_moving",t,".pdf"),width=12,height=16)
plotAllTrajVar(moving_theta[moving_theta$vt == rates[t],],m=0.1,E=0,obs="dist")#,ylim=c(0,.005))
#dev.off()
}
Check and compare the final variance :
for(t in seq_along(rates)){
#pdf(paste0("equilibrium_var_x_moving",t,".pdf"),width=12,height=16)
plotAllVariableSummaries(moving_theta[moving_theta$vt == rates[t],],E=0,estimate=eq2833b)
#dev.off()
}